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Abstract 

We present results from a study of the nonlinear intermodal coupling between different flexural 
vibrational modes of a single high-stress, doubly-clamped silicon nitride nanomechanical beam. The 
measurements were carried out at 100 mK and the beam was actuated using the magnetomotive 
technique. We observed the nonlinear behavior of the modes individually and also measured the 
coupling between them by driving the beam at multiple frequencies. We demonstrate that the 
different modes of the resonator are coupled to each other by the displacement induced tension in 
the beam, which also leads to the well known Duffing nonlinearity in doubly-clamped beams. 
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I. INTRODUCTION 



The nonlinear dynamics of nanoelectromechanical systems (NEMS) has attracted con- 



siderable interest over recent years. 



lHl2j Nonlinear behavior is of practical importance in 



NEMS as it can limit their use as sensors operating in the linear regime, [4l 1 1 31 ] though it 
also opens up the possibility of devices that are designed to exploit it.j^, Q However, the 
very high operating frequencies of NEMS devices (which typically have resonant frequen- 
cies in the MHz range) make them ideal for investigating fundamental aspects of nonlinear 



dynamics. 



12 1 Because NEMS devices are engineered systems and their nonlinear- 

n 

ities can be extensively tuned, |6[ it should be possible to use them to study complex non- 



linear phenomena 



15| in which a large number of different mechanical modes interact 



Furthermore, when cooled to very low temperatures, the nonlinear behavior in nanomechan- 
ical devices would provide important signatures of the transition from classical to quantum 
regimes. Il7f 

n 

Nonlinearity can arise in a wide variety of different ways in NEMS devices. 2 External 
potentials generated electrostatically by fixed electrodes can be used to induce (or suppress) 
nonlinear behavior in NEMS.|(| However, nonlinearities can be intrinsic to a given device or 
material; for example when a doubly clamped beam is excited its displacement leads to a 
stretching of the beam which generates a Duffing type nonlinearity in a given mode as well 
as generating inter-mode couplings. 10J, The damping of nanomechanical devices can also 
be quite nonlinear, something which appears to be particularly important for carbon based 



devices. 



Experiments have begun to explore the coupled nonlinear dynamics of mechanical modes 



n 



either in two adjacent resonators |9j or occuring as different harmonics within a single 



beam. 
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ll| In a recent work Westra et al. 10] explored the nonlinear coupling between the 



first and third harmonics of a micromechanical doubly clamped beam with mode frequencies 
in the kHz range. Although the Q-factor of the modes was rather low (~ 10 2 ) they were 
able to see clear evidence of nonlinear coupling between the modes induced by the stretching 
nonlinearity in their beam. 

In this paper, we present a systematic study of the nonlinear response of three of the flex- 
ural modes of a nanomechanical beam resonator, together with the associated inter-mode 
couplings. We report the results of a series of measurements taken at 100 mK on the funda- 
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mental, third and fifth harmonics (all of which are in the MHz range) of a doubly-clamped 
beam fabricated from high-stress (~ 1 GPa) silicon nitride. In a first set of measurements, 
we excite one mode at a time and find that they display a Duffing like behavior which can 
be understood quite naturally in terms of the nonlinearity arising from the stretching of 
the beam. We then measured the response when two modes are driven at the same time, 
with one mode acting as a probe of the displacement of the other. The high tensile stress 

n 

and the low temperature mean that the Q-factor of the modes is high enough jl9| to allow 
us to observe the higher mode resonance signals as well as the frequency shifts induced 
by nonlinear mode-mode couplings which are much larger than the linewidth. Using the 
single-mode results to calibrate our measuring scheme we are able to make a detailed quan- 
titative comparison (without the need for any fitting parameters) between our results and 
the predictions of a simple model of the stretching nonlinearity. 

This paper is organized as follows. Section II introduces a simple theoretical model of 
the stretching nonlinearity in beam resonators. The experimental set-up is outlined in Sec. 
Ill and in Sec. IV we describe our results. We give our conclusions in Sec. V. 



II. THEORETICAL BACKGROUND 



Here we review the theoretical description of a beam under tension, 20[ starting from 
Euler- Bernoulli equation extended to include the geometric nonlinearity arising from the 
stretching of the beam that accompanies its flexure. 0, [3] For a single mode, this non- 
linearity leads naturally to the Duffing equation where a term proportional to the cube 
of the displacement is added to the usual harmonic equation of motion.j2l However, the 
non-linearity also generates important couplings between different modes. [10] 

We start by considering a doubly-clamped beam of length L, lying along the x-axis, with 
cross-sectional area A = wh, where w is the width and h the thickness. The beam is assumed 
to be under intrinsic tension, T , and we include, to lowest order, the nonlinear tension arising 
from the stretching of the beam. Q The equation of motion for the displacement, y, is given 
by 



pAy + V y + EI y y h 



Tn 



EA 
~2L 



(y) dx 



y 



(i) 



where p is the density, r\ characterizes the damping, I y = hw 3 / 12 is the moment of inertia, 
Fl is the force per unit length exerted on the beam and E the Young's modulus. The 
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non-linearity arises from the second term within the square brackets. 

The mode frequencies and mode functions for the corresponding linear problem are 
known 20| and provide a convenient starting point for describing the non-linear regime. 
If the system is driven harmonically, Fl = F n cos(u n t), at a frequency close to one of 
the linear mode frequencies w n ,o> we can then approximate the displacement of the beam 
as y(x,t) = wg n (x)u n (t) where g n (x) is the n-th mode function, normalized so that 
Jq 1 g 2 (x)dx = 1 where x = x/L. 

The dimensionless function u n (t) obeys the equation of motion,^ 



+ -TT 7- + W n A + ^nn U n = In COs(w n t), 



(2) 



dt 2 ' Q n dt 

where A = (Ew 2 )/(2pL 4 ), f n = F n £ n /(hw 2 p), with £ n = J dxg n (x) (the mode parameter), 
and the Q-factor of the n-th mode is given by Q n = u n ^{Ap/rj). The strength of the 
nonlinearity is controlled by the parameter I nn which is defined by an integral of the form, 
Iij = dxgi(x)g'Ax) (with i and j labeling any two modes of the beam). 

Since the equation of motion [Eq. [2] is that of the familiar Duffing oscillator, we proceed 
to analyse it using standard methods. |2lj] Substituting a solution of the form, 



a n cos(u n t) + b n sin(u;„t), 



into Eq. [21 we obtain an equation for the amplitude of the motion, r n = (a^ + fe^) 1 ^ 2 , 



fn 



U n,0 + ^4n r n ~ UJ n) + 



Qn 



-1/2 



(3) 



(4) 



For a large enough quality factor, Q n , and for relatively small amplitudes, r n , the main effect 
of the nonlinearity is to shift the frequency of peak response to 



(5) 



In the weakly nonlinear regime, the frequency of the resonator shifts from its natural (un 



driven) value quadratically with increasing amplitude (frequency pulling). 2l|, [22 1 



We now consider a situation where the drive contains two harmonic components, 

F L = F n cos(u n t) + F m cos(u m t), (6) 

with the two drives frequencies u n and u m chosen to be close to the resonant frequencies of 
two different modes, uq u and uo m , respectively. In this case we assume the beam displace- 
ment can be approximated as, y(x,t) = w[u n (t)g n (x) + u m (t)g m (x)]. 



Following the same approach as for the one mode case, assuming u n and u m each oscillate 
at a single frequency (i.e. they take the form of Eq. [3]), we obtain a modified expression for 
the amplitude of mode n, 

-1/2 



fn fn 

where 



/ \ 2 

I 2 2\2 . / W n U„,o 

K,2 - U n ) + [—q- 



(7) 



u n ,2 = w n ,o\l 1 + ^A/2 n r2 + A Q J nn I mm + i^J r 2 m . (8) 

Thus we see that the frequency of a given mode depends not just on the amplitude of its own 
motion, but the amplitude of the other mode that is excited. In particular, the frequency 
shift of one mode initially grows quadratically with the amplitude of the other one. 

III. EXPERIMENTAL SET-UP 

Experiments were performed on a silicon nitride beam with length 25.5 /iin, width 170 nm 
and thickness 170 nm. Metal layers consisting of 3 nm of Ti and 40 nm of Au were added by 
thermal deposition to form a wire on top of the beam through which a drive current could 



be applied. The devices were fabricated by dry-etching in a multi-stage process 26j] from 



wafers composed of a 190 fim thick silicon wafer with 170 nm of silicon nitride on both sides; 



the nitride layer has a built-in tensile stress of about 1070 MPa at room temperature. 23] 

A schematic diagram of the experimental set-up is shown in Fig. [Th- The resonator 
was placed inside a 3 He/ 4 He dilution refrigerator and measurements were performed at 



100 mK. The motion of the beam was monitored using the the magnetomotive method. [241] 
A magnetic field of B = 3 T was applied and the drive was applied by passing an alternating 
current through the sample. The response to the resulting Lorentz force was detected by 
measuring a voltage (V a ), which is related to the input voltage at the rf pre-amplifier (Vb), 
as the frequency of the drive signal was tuned through the mechanical resonance. The beam 
was driven at two different frequencies by using an additional signal generator as shown in 
Fig. [1^. The measurements were carried out in transmission mode, see Fig. [lb, capturing a 
dip in the conductance of the beam in the spectral domain. 

The motion of the resonator in the magnetic field generates an emf,[24] which itself is 
related to the mechanical motion. Provided that the resonances are well-resolved, \uo m — 
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FIG. 1: (a) Schematic diagram of the measurement set-up. The different modes of the beam were 
driven with two separate signal generators. A voltage signal proportional to the emf generated by 
the beam was detected at the input of the network analyzer after passing through a pre-amplifier. 
(b) Circuit model for transmission measurements, with external circuit impedances Z e and Zq, 
sample impedance Z m , the drive voltage V m and the measured voltage Vq. 

w on| ^ tt 21 + ttN for m 7^ n, the motion at a given frequency can be attributed to a single 
mode and 

V cmi = w£ n LB^. (9) 

We can relate the final measured voltage, V an to Kmf> an d hence to the amplitude of a given 
mode, r n , by applying Kirchoff 's law to the circuit (Fig. [T|d) and solving the corresponding 
equations. Thus we find, 

V an = GV = (V in - V cmi ) , (10) 

z + z e 
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FIG. 2: Resonant response of mode 1 (7.5 MHz). The curves from left to right correspond to a 
series of increasing drive strengths. The normalized response consists of the measured voltage, V a i, 
scaled by the maximum value measured in each case. 



where G is a constant of proportionality that quantifies the overall signal gain from the 
sample to the network analyzer, taking into account the losses in the transmission cables 
as well as the gain provided by the amplifiers. The gain in the transmission line will in 
fact vary slightly with frequency (hence the gain will not be the same for different modes) 
and local temperature variations inside the fridge. Mechanical resonance leads to a dip in 
the voltage V an , hence (assuming the circuit impedances are real) we can characterize the 
mechanical response of a particular mode by, 

Z G 



Sn 



z + z, 
z 



-Kn - V„, 



GBLw^ n uj n r n , 



(11) 



which measures the size of the dip. Expanding Eq. [5] to lowest order in r n and using Eq. [TT] 
leads to the explicit expression: 



where a r 



3I^ n \Ld„fiDn 

S,(BLu n iw) 2 



and D n 



( Z + Ze 1 

\ Z G 



(12) 



7 



IV. RESULTS 



We were able to detect the first three odd harmonics of the beam which had frequencies 
7.5, 22.85 and 39.28 MHz respectively (the even modes do not give rise to a signal which is 
directly measurable using the magnetomotive method). Measurements on each individual 
mode were performed first, capturing the response of the resonator as the driving frequency 
was swept through the mechanical resonant frequency. The basic properties of the three 
modes are summarized in Table I. The frequencies are those measured at the lowest drive 
levels used. The Q-factors are the measured values at a field of 3 T. Although the values are 
still rather high, the relatively strong magnetic field used means that they are lower than 
the intrinsic Q-factors of the modes [24J (much lower in the case of the n — 1 mode 26j whose 
intrinsic Q-factor is 1.8 x 10 6 at 100 mK). 

The nonlinear parameters, I nn , and the mode parameters, £ n , are calculated numerically 
using a value of the intrinsic tension which we estimate to be To = 1020 MPa. This is slightly 
lower than the room temperature value because of the different thermal contractions of the 



silicon and silicon nitride layers in the wafer. 



26] The parameters, I nn , and hence the strength 



of the nonlinearity in a given mode, grow rapidly in size with the mode number n. The cross 
terms, I nm , which involve two different modes, also grow with the mode number, but are 
much smaller for our device: /13 = —1.9, /15 = —3.0, and J35 = —8.8. In calculating the 
parameter A = (Ew 2 ) / (2pL 4 ) we use the Young's modulus of silicon nitride {27] E S iN = 211 
GPa and neglect the contribution from the gold as it is under much less tension and its 
Young's modulus, E Au = 78 GPa, is smaller. We include the effect of the gold layer on the 
mass by using an appropriate average value for the density. Using these parameters we obtain 
theoretical estimates of the three mode frequencies of 7.8, 24.2 and 42.5 MHz respectively 
which are all close to the measured values (given in Table 1). The discrepancies between 
estimated and measured frequencies arise from the uncertainties in the device dimensions 
(L ~ ±1%, w ~ ±5% and t Au ~ ±10%). 

Data were collected for a range of drive amplitudes by systematically increasing the drive 
signal. As an example, the spectral response of the fundamental mode is shown in Fig. [2] as 
a function of the corresponding drive amplitude. As expected, the curves are symmetric for 
the smallest drives, but become increasingly asymmetric as the beam is driven harder. At the 
largest amplitudes, the device enters a strongly nonlinear regime marked by sharp changes 
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FIG. 3: Shift in the frequency, Afn > as a function of the voltage Vs n for (a) n = 1, (b) n = 3, 
and (c) n = 5. The points are the data and the lines are quadratic fits. The corresponding values 
of the parameters D n obtained from the fit, together with the associated uncertainty, are given in 
each case. 



n 


a;o„/2vr (MHz) 




Inn 


Qn 


1 


7.50 


0.88 


10.4 


9.0 x 10 3 


3 


22.85 


0.30 


93.3 


1.7 x 10 5 


5 


39.28 


0.19 


257.9 


4.1 x 10 5 



TABLE I: Properties of the three measured modes of the silicon nitride beam. 



seen in the measured signal on the high frequency side of the resonance. 28j Frequency 
pulling is also clearly visible, with the peak frequency of each mode shifting upwards as the 
drive is increased. Very similar behavior is seen for modes 3 and 5. 

We extracted the frequencies of peak response, ffi = oj n ^/2n and the corresponding 
voltages Vsn from the measured resonance curves (Fig. [2]). The shift in peak frequency, 
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FIG. 4: Response curves of the third flexural mode for different measured voltages, Vsi, at the 
frequency of the first mode, V\ = 15.2 mV, V2 = 29.3 mV and V3 = 48 mV, corresponding to three 
successively increasing drive amplitudes applied at a fixed frequency of 7.499 MHz. 

A/n = ( w n,i — ^n,o)/27r, versus peak measured voltage Vg m for each mode is shown in Fig. 
EJ Theory predicts (see Eq. ITS!) that the frequency shift should increase quadratically with 
the voltage and this is what is found. All the parameters required for the analysis are known 
(with stated uncertainty) except for the D n . We obtain a value of the quadratic coefficient, 
a n , for each mode by fitting each set of measured data to a quadratic function with a fitting 
error of less than 3%, from which we then calculate the D n . The fits are shown as lines in 
Fig. |3j The value of D n varies between the modes by about a factor of two, due to the 
intrinsic frequency dependence of the circuit impedances and amplifier gain. The error in 
D n is about 13% and follows directly from the error in sample dimensions. 

Next we probed the interactions between the different modes of the beam by exciting one 
mode weakly and measuring its response as successively stronger drive amplitudes (of a fixed 
frequency) were applied to a second mode. The shift in the peak response frequency of the 
weakly driven mode, A/^ = (uj n ^ — w ni o)/2vr, was measured together with the voltage at 
the frequency of the second mode. Response curves for mode 3, measured for three different 
levels of drive applied to mode 1, is shown in Fig. HI As expected, the resonant frequency 
of mode 3 increases with the measured voltage from the first mode. Looking at the depths 



of the curves in Fig. HI we see that increasing t 
the amplitude of the response of mode 3 slightly. 



le amplitude of the mode 1 drive reduces 



1 Of ] A reduction in the peak amplitude is a 
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FIG. 5: Frequency shifts of (a) the fifth mode, A/g , and (b) the third mode A/3 , as a function 
of the measured response of the first mode, Vsi, to a sequence of progressively stronger drives. In 
each case the points are data, and the shaded area is the band of values consistent with theoretical 
predictions when uncertainties in the parameters are accounted for. 

natural consequence of an increase in the frequency u n> 2 [see Eq. [7] , but we also find a very 
slight degradation in the Q-factor of the third mode accompanies the increasing drive at the 
fundamental frequency. 

The shift in the peak response frequencies of the third and fifth modes as a function of 
the voltage measured at the frequency of mode 1, V51, are shown in Fig. |5j Similarly, Fig. 
[6] shows the effect of varying the amplitude of the third mode on the resonant frequency of 
the first and fifth modes. Again the theory predicts (through Eqs. M and (HI) a quadratic 
dependence of the frequency shifts on the measured voltages. However, because the values of 
the parameters a n have been determined from fits to the single-mode data, the comparison 
with theory now involves no free parameters. We use the a n values extracted from the 
single mode data, and the parameters I nm , together with the associated uncertainties, to 
plot shaded bands in Figs. and [H] showing the regions which are fully consistent with the 
theory. We note that the I nm values are quite insensitive to the changes in beam dimensions, 
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FIG. 6: Frequency shifts of (a) the fifth mode, A/g , and (b) the first mode A/{ , as a function 
of the measured response of the third mode, V53, to a sequence of progressively stronger drives. In 
each case the points are data, and the shaded area is the band of values consistent with theoretical 
predictions when uncertainties in the parameters are accounted for. 

so the predictions of the theory are fairly precise. In each case it is clear that the dependence 
of the frequency shifts on the voltages is well described by a quadratic law, and that there 
is very good agreement between the theoretical predictions and the measurements. 

Having verified that the bending nonlinearity provides a quantitatively accurate descrip- 
tion of the modal couplings in the beam, we can now obtain the amplitude of the motion in 
one mode by measuring the frequency shifts in a second (weakly driven) mode. The theo- 
retical expression, Eq. [HJ and the mode function g n (x) allow us to convert from a measured 
frequency shift to a physical amplitude at a given point along the beam. An example is 
shown in Fig. [7] in which the amplitude of mode 1 at the antinode x = L/2, obtained by 
measuring a frequency shift in mode 3, is shown as a function of the drive frequency, 
An amplitude of lnm in mode 1 translates into a frequency shift in mode 3 of 7.11Hz, a 



factor of ~ 40 larger than that observed in Ref . QJJ. From Fig. [7J we see that an amplitude 



of 25 nm (~ 10% of the width of the beam) is already well within the non-linear regime. 
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FIG. 7: Amplitude of the motion of mode 1 (at the antinode), A±, as a function of the drive fre- 
quency wi/2-7r obtained using the frequency shift in mode 3. Each point is generated by measuring 

(2) 

a full frequency response curve for mode 3 and the obtaining the resulting A/3 which is then 
converted into an amplitude. 

V. CONCLUSIONS 

We have measured the properties of three of the flexural modes of a highly stressed, 
doubly-clamped silicon nitride nanomechanical beam. The size of the sample, the Q-factors, 
the temperature and the degree of non-linearity are all orders of magnitude different than 



reported by Westra et al. 10[ for an unstressed beam at room temperature. As a first step 
we investigated the frequency pulling that accompanies increases in the drive applied to each 
of the modes in turn. For all three modes the frequency was found to grow quadratically 
with the amplitude of the motion, consistent with the Duffing-type behavior predicted for 
a nonlinearity governed by the stretching of the beam on deflection. 

We then examined the behavior of the system when one mode is driven weakly and the 
amplitude of the drive applied to a second mode is increased steadily. We again found 
that the frequency of the weakly driven mode increased quadratically, this time with the 
amplitude of the second mode. Using the results from the single-mode experiments to 
calibrate our measurement set-up we were able to make a comparison without any free 
parameters to theoretical predictions. The good agreement we find using four different pairs 
of modes allows us to conclude that the nonlinear dynamics we observe does indeed arise 
from the stretching of the beam on deflection. 
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In conclusion, our data confirm that the bending nonlinearity provides a good quantitative 
description of the mode couplings in our nanomechanical device. We find that beam starts 
to behave nonlinearly for amplitudes that are a small fraction of its width, which means that 
the system has a high intrinsic nonlinearity due to its size and the built-in stress in the nitride 
layer. Furthermore, our measurements also allows us to calculate the physical displacement 
of one mode by measuring a frequency shift in a second mode. Studying the inter-modal 
couplings in nanomechanical beams is of great importance for understanding the dynamics 
of such small systems, and can be of great use to device applications or fundamental studies 
which require a combination of high frequencies, large Q-factors and sensitivity. 
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